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ABSTRACT 

We obtain predictions for the properties of cold dark matter annihilation radiation using 
high resolution hydrodynamic zoom-in cosmological simulations of Milky Way-like galax¬ 
ies (apostle project) carried out as part of the “Evolution and Assembly of GaLaxies and 
their Environments” (EAGLE) programme. Galactic halos in the simulation have signihcantly 
different properties from those assumed in the “standard halo model” often used in dark mat¬ 
ter detection studies. The formation of the galaxy causes a contraction of the dark matter halo, 
whose density profile develops a steeper slope than the Navarro-Erenk-White (NEW) profile 
between r « 1.5 kpc and r « 10 kpc. At smaller radii, r < 1.5 kpc, the halos develop a 
flatter than NEW slope. This unexpected feature may be specihc to our particular choice of 
subgrid physics model but nevertheless the dark matter density prohles agree within 30% as 
the mass resolution is increased by a factor 150. The inner regions of the halos are almost 
perfectly spherical (axis ratios b/a > 0.97 within r = 1 kpc) and there is no offset larger 
than 45 pc between the centre of the stellar distribution and the centre of the dark halo. The 
morphology of the predicted dark matter annihilation radiation signal is in broad agreement 
with 7 -ray observations at large Galactic latitudes (b > 3°). At smaller angles, the inferred 
signal in one of our four galaxies is similar to that which is observed but it is signihcantly 
weaker in the other three. 

Key words: cosmology: theory, dark matter, gamma-rays: galaxies. Galaxy: centre, methods: 
numerical 


1 INTRODUCTION 

Uncovering the unknown nature of the dark matter is one of the 
greatest challenges of modem cosmology and particle physics. 
Since the original suggestion that the dark matter might consist of 
massive, cold, weakly interactive, neutral particles, a large body of 
empirical evidence has consolidated this hypothesis which, how¬ 
ever, can only be confirmed by the detection of the particles them¬ 
selves. Among the possible candidates, supersymmetric particles 
(see Jungman et al. 1996, for a review) such as neutralinos are at¬ 
tractive options that current particle accelerator experiments might 
detect. 

An interesting property of many particle candidates for cold 
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dark matter is that they annihilate into standard model particles, in¬ 
cluding photons. This opens up the exciting possibility of attempt¬ 
ing to detect such photons from space. The requirement that weakly 
interacting particles provide the measured dark matter density in 
the Universe today suggests a plausible particle mass range of or¬ 
der = 10 — 1000 GeV, leading to the emission of 7 -ray pho¬ 
tons in or below that energy range when two dark matter particles 
annihilate (see review by Bertone et al. 2005). 

The Large Area Telescope aboard the Gamma-Ray Space 
Telescope (Fermi) (Gehrels & Michelson 1999) has, over the last 
few years, produced the most detailed maps of the 7 -ray sky, cov¬ 
ering a large energy range (20 MeV — 500 GeV) with a resolu¬ 
tion of a few arcmin at the highest energy end of the spectrum. 
Analysing the Fermi data around the Galactic Centre (GC), a num¬ 
ber of authors (Goodenough & Hooper 2009; Hooper & Good- 
enough 2011; Hooper & Linden 2011; Abazajian & Kaplinghat 
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2012; Gordon & Macias 2013; Hooper & Slatyer 2013; Macias & 
Gordon 2014; Daylan et al. 2014; Abazajian et al. 2014; Galore 
et al. 2015b) have claimed the detection of extended diffuse excess 
emission above the other known astrophysical sources. This excess 
emission, peaking at i? ~ 2 GeV, was found to be broadly consis¬ 
tent with the expected signal from dark matter annihilation. In par¬ 
ticular, the flux decreases with distance from the GC as with 
slope, r = 2.2 — 2.4, only slightly steeper than the asymptotic in¬ 
ner slope, r = 2, of flux originating from the NFW (Navarro et al. 
1996b, 1997) density profiles found in dark matter only simulations 
of halos. 

However, potential systematic effects in the analysis of the 
7 -ray data could be introduced by incorrect point source subtrac¬ 
tion or the modelling of the diffuse backgrounds (Ackermann et al. 
2012). Alongside the dark matter interpretation of the Galactic 
Centre excess, astrophysical explanations have been proposed. For 
example, a population of as yet unresolved millisecond pulsars 
(MSP) (Abazajian 2011; Gordon & Macias 2013; Yuan & Zhang 
2014; Cholis et al. 2015b) or young pulsars (O’Leary et al. 2015) 
associated with the bulge^ of the Milky Way. In contrast to the con¬ 
clusions of Hooper et al. (2013) and Cholis et al. (2015b), Bar¬ 
tels et al. (2015) and Lee et al. (2015) have argued that the ex¬ 
cess could be due to such a pulsar population. Alternatively, Linden 
et al. (2012); Carlson & Profumo (2014); Petrovic et al. (2014) and 
Cholis et al. (2015a) have suggested that the excess could origi¬ 
nate from the injection of very high energy cosmic rays during past 
activity in the Galactic Centre. 

Besides the spectral shape, another property that can help dis¬ 
tinguish the potential sources of 7 -rays contributing to the excess 
is the morphology of the signal. A dark matter origin would re¬ 
quire the excess to extend over tens of degrees on the sky (Springel 
et al. 2008b; Setpico & Zaharijas 2008; Nezri et al. 2012). An ex¬ 
cess with the same spectral shape extending over a large angu¬ 
lar range, with emissivity decreasing with distance from the GC, 
would strengthen the interpretation of the excess as originating 
from dark matter annihilation. Using multiple regions between 2° 
and 20° from the GC and a large range of Galactic diffuse emission 
(GDE) models, Calore et al. (2015b) found that the excess emission 
is consistent with a dark matter particle of mass GeV 

annihilating into a bb quark pair then producing photons and is dis¬ 
tributed following a generalised NFW profile (gNFW, see equation 
1 below) with a slope 7 = 1.26 ±0.15.A similar spatial distribu¬ 
tion was found by Daylan et al. (2014) who suggested a slope for 
the inner profile in the range 7 = 1.1 — 1.3. With the increasing 
precision of these measurements and of the foreground modelling, 
it has become crucial to refine the theoretical models for the distri¬ 
bution of dark matter at the centre of galaxies in a ACDM context. 
Characterising the dark matter profile slope and sphericity as well 
as investigating the potential offset between the dark matter and the 
GC are all important tasks for theorists. 

Work based on dark matter only simulations has shown that 
the dark matter is distributed following an NFW density profile 
with a scale length, r^, which varies with halo mass (e.g. Navarro 
et al. 1997; Neto et al. 2007; Duffy et al. 2008; Dutton & Maccio 
2014). Higher-resolution simulations (Springel et al. 2008a) have 
shown that the very inner parts of dark matter profiles might be 
slightly shallower than the asymptotic NFW slope of 7 = 1 
(Navarro et al. 2010). Similarly, predictions for the signal com- 
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ing from subhaloes have also been made using these simulations 
(Kuhlen et al. 2008; Springel et al. 2008b), effectively proposing 
a test of the cold dark matter paradigm. At the other end of the 
halo mass range, Gao et al. (2012) argued that nearby rich clusters 
provide a signal with a higher signal-to-background ratio than the 
Milky Way’s satellites. Thus far observational measurements have 
proved inconclusive and the only claimed detection comes from the 
centre of our own Milky Way, where precise predictions from dark 
matter simulations have been made (Springel et al. 2008b). 

However, these studies all ignored the effects of the forming 
galaxy on the structure of the dark matter halo. Mechanisms such 
as dark matter contraction (e.g. Barnes & White 1984; Blumen- 
thal et al. 1986; Gnedin et al. 2004) can drag dark matter towards 
the centre, steepening the profile. Conversely, perturbations to the 
potential, due for instance to feedback from stars or supermassive 
black holes or the formation of a bar, can lead to a flattening of the 
very central regions (e.g. Navarro et al. 1996a; Weinberg & Katz 
2002; Mashchenko et al. 2006; Pontzen & Govemato 2012). The 
correct balance between these mechanisms can only be understood 
by performing high-resolution hydrodynamic simulations of Milky 
Way-like galaxies using a physical model validated by comparison 
with a wide range of other observables. 

In this study we use two “zoom” simulations of Local Group 
environments (the APOSTLE project: Sawala et al. 2015; Fattahi 
et al. 2015) performed within the framework of the “Evolution 
and Assembly of GaLaxies and their Environments” (EAGLE) suite 
(Schaye et al. 2015; Crain et al. 2015). These simulations have been 
shown to reproduce a large number of observables of the galaxy 
population at low and high redshifts, as well as the satellite galaxy 
luminosity functions of the Milky Way and Andromeda galaxies 
with unprecedented accuracy. Schaller et al. (2015a) showed that 
the EAGLE simulations produce galaxies with rotation curves that 
are in unprecedented agreement with observations of field galaxies, 
suggesting that the matter distribution in the simulated galaxies is 
realistic and that the main that effects of baryons on halos are accu¬ 
rately captured by the simulations. Note, however, that Calore et al. 
(2015a) showed that the goodness of fit of the simulated data to 
the observed Milky Way rotation curve is lower in the highest res¬ 
olution zoomed-in simulations. The EAGLE simulations therefore 
provide an excellent test-bed for the interpretation and analysis of 
the Fermi excess. 

This paper is structured as follows. In Section 2 we introduce 
the simulation setup used. We investigate the dark matter density 
profiles of our halos in Section 3 and analyse the dependence of the 
annihilation signal on the angle with respect to the Galactic Centre 
in Section 4. We summarise our findings and conclude in Section 
5. 

Throughout this paper, we assume a WMAP7 flat ACDM cos¬ 
mology (Komatsu et al. 2011) {h = 0.704, Qh = 0.0455, 11^ = 
0.272 and ag = 0.81), express all quantities without h factors and 
assume a distance from the GC to the Sun of tq = 8.5 kpc. 


2 THE SIMULATIONS 

The simulations used in this study are based on the EAGLE simu¬ 
lation code (Schaye et al. 2015; Crain et al. 2015) . We summarize 
here the parts of model relevant to our discussion. 
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2.1 Simulation code and subgrid models 


The EAGLE code is based on a substantially modified version of 
the GADGET code, last described by Springel (2005). The mod¬ 
ifications include the use of a state-of-the-art implementation of 
Smoothed Particle Hydrodynamics (SPH), called ANARCHY (Dalla 
Vecchia (in prep.), see also Schaller et al. 2015c), based on the 
pressure-entropy formulation of SPH (Hopkins 2013). Gravita¬ 
tional interactions are computed using a Tree-PM scheme. 


The cooling of gas and its interaction with the background 
radiation is implemented following the recipe of Wiersma et al. 
(2009a) who tabulated photoheating and cooling rates element-by- 
element (for the 11 most important elements) in the presence of the 
UV and X-ray backgrounds inferred by Haardt & Madau (2001). 
To prevent artificial fragmentation, a pressure floor in the form of 
an effective equation of state, Peos oc , designed to mimic 
the mixture of phases in the interstellar medium (ISM) (Schaye & 
Dalla Vecchia 2008), is applied to the cold and dense gas. Star for¬ 
mation is implemented using a pressure-dependent prescription that 
by construction reproduces the observed Kennicutt-Schmidt star 
formation law (Schaye & Dalla Vecchia 2008) and uses a density 
threshold that captures the metallicity dependence of the transition 
from the warm, atomic to the cold, molecular gas phase (Schaye 
2004). Star particles are treated as single stellar populations with 
a Chabrier (2003) initial mass function (IMF) evolving along the 
tracks advocated by Portinari et al. (1998). Metals from supernovae 
(SNe) and AGB stars are injected into the ISM following the model 
of Wiersma et al. (2009b) and stellar feedback is implemented via 
the stochastic injection of thermal energy into the gas neighbour¬ 
ing newly-formed star particles as described by Dalla Vecchia & 
Schaye (2012). Galactic winds hence form naturally without having 
to impose a direction, velocity or mass loading factor. The amount 
of energy injected into the ISM per feedback event is dependent 
on the local gas metallicity and density in an attempt to take into 
account the unresolved structure of the ISM (Schaye et al. 2015; 
Crain et al. 2015). Supermassive black hole seeds are injected in 
halos above 10 ^°/i~^Mq and grow through mergers and the ac¬ 
cretion of low angular momentum gas (Rosas-Guevara et al. 2015; 
Schaye et al. 2015). AGN feedback is modelled by stochastically 
injecting thermal energy into the gas directly surrounding the black 
hole (Booth & Schaye 2009; Dalla Vecchia & Schaye 2012). Halos 
are identified using the FOF algorithm (Davis et al. 1985) and sub¬ 
structures within them are identified in post-processing using the 
SUBFIND code (Dolag et al. 2009). 


The subgrid model was calibrated (by adjusting the efficiency 
of stellar feedback and the accretion rate onto black holes) so as to 
reproduce the present day stellar mass function and galaxy sizes, as 
well as the relation between galaxy stellar masses and supermas¬ 
sive black hole masses (Schaye et al. 2015; Crain et al. 2015). As 
argued by Schaye et al. (2015) (Sec. 2), numerical convergence in 
the traditional sense (strong convergence) cannot be achieved when 
new physical processes are resolved at each resolution. In this case, 
the free parameters of the model should be re-calibrated to match 
the same pre-defined set of observables (weak convergence). This 
can be done in cosmological simulations of representative periodic 
volumes (Crain et al. 2015), but it is much more difficult to achieve 
for “zoom-in” simulations of a few objects. In this case, even weak 
convergence is difficult to establish (See Sec. 3.3). 
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2.2 Selection of Milky Way halos 

The two volumes used in this work are zoom resimulations of re¬ 
gions extracted from a dark matter only simulation of 100 ^ Mpc® 
with 1620® particles. The halos were selected to match the 
observed dynamical constraints of the Local Group (APOSTLE 
project: Sawala et al. 2015; Fattahi et al. 2015). Each volume con¬ 
tains a pair of halos in the mass range M 200 = 5 x 10 ®®Mq to 
M '200 = 2.5 X 10®®Mq that will host analogues of the Milky 
Way and M31. We use the two halos in volumes AP-1 and AP-4 
of the APOSTLE suite (see Table 2 of Fattahi et al. 2015, where 
other relevant data are listed). The high-resolution region encloses 
a sphere larger than 2.5 Mpc around the centre of mass of the two 
halos z — 0. The dark matter particle mass in the zoom re¬ 
gions was set to 5 X 10"^Mq, whilst the primordial gas particle 
mass was set to 1 x IO'^Mq. The initial conditions were generated 
from ACDM power spectra using 2"® order Lagrangian perturba¬ 
tion theory (Jenkins 2010) and linear phases were taken from the 
Gaussian white noise field PANPHASIA (Jenkins 2013). The gravi¬ 
tational softening length was set to e = 134 pc (Plummer equiv¬ 
alent) at 2 < 2.8 and was kept fixed in comoving units at higher 
redshifts. Simulations without baryonic components were run with 
the exact same setup and are labelled as DMO in what follows. 


3 DARK MATTER DISTRIBUTION IN THE CENTRE OE 
THE HALOS 

In this section, we analyse the dark matter distribution of the haloes 
of simulated Milky Way galaxies. We consider both central galaxies 
in each of the two simulation volumes as Milky Way-like galaxies. 


3.1 Profiles without haryon effects 

The analysis of the GC excess is often performed using an assumed 
analytic density profile shape for the dark matter. This profile is a 
generalisation of the NFW profile for which the asymptotic inner 
slope is a free parameter 7 : 


Pdm(t) 


Ps _ ^ 

{r/rf^ (1 + 


( 1 ) 


The NFW profile is recovered for 7 = 1 . This generalized form 
of the density profile is not supported by numerical simulations of 
collisionless dark matter (e.g. Navarro et al. 2010) but is a useful 
way to parametrise the deviation from the NFW shape in the very 
centres of halos as a result of baryonic effects. As the measurements 
of the GC excess only span a range of a few kiloparsecs, the value 
of the scale radius Va cannot be constrained observationally and is 
typically fixed to Ts = 20 kpc, in broad agreement with simulation 
results for Milky Way-like halos (e.g. Neto et al. 2007; Dutton & 
Maccio 2014). The normalisation of the profile, ps, is degenerate 
with other particle physics parameters (see Section 4) and is usually 
fixed by requesting that the density of dark matter at the location of 
the Sun^ is Pdm(tq) = 0.4 GeV • cm“®, in agreement with local 
dynamical constraints (Catena & Ullio 2010; locco et al. 201 1). 

In order to quantify the effects of baryons on the dark matter 
distribution, it is worth first considering the profiles extracted from 


^ Note that for simplicity we use units convenient for particle physics ap¬ 
plications. Units more friendly to astronomers are recovered using the con¬ 


version 1 Mq • kpc ® = 3.795 X 10 ® GeV ■ cm 


-3 
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Figure 1. The dark matter density profiles of the four halos in the simula¬ 
tions without baryons (yellow solid lines). Thinner lines are used at radii 
smaller than the convergence radius of the simulation. The vertical dotted 
line indicates the simulation’s gravitational softening length. The green dot- 
dashed and blue dashed lines comespond to an NFW and gNFW with 7 = 
1.26 profiles respectively, both normalised to p(rQ) = 0.4 GeV ■ cm“® 
and with a scale-length Ts = 20 kpc. As expected, the simulated profiles 
display a shape similar to the plotted NFW model, but with a lower normal¬ 
isation than the standard halos. 


the simulations without baryonic physics. In Fig. 1, we show the 
dark matter density profiles of our four halos. Thick lines are used 
at radii greater than the resolution limit (Rpos ~ 350 — 450 pc 
depending on the halo) set by the criterion of Power et al. (2003) 
and thin lines are used at smaller radii. The softening length is in¬ 
dicated by the vertical dotted line. The green dot-dashed and blue 
dashed lines correspond to NFW and gNFW with 7 = 1.26 (the 
best-fitting value of Galore et al. (2015b) to the excess) profiles re¬ 
spectively, both normalised, as discussed above, to pDM(f 0 ) = 
0.4 GeV • cm“® and with a scalelength rs = 20 kpc. As expected, 
the profiles are in good agreement with the NFW model albeit with 
a lower normalisation. The best-fitting NFW profile parameters to 
our halos are given in Table 1. The usual choice of = 20 kpc 
is in good agreement with our simulated halos but the normalisa¬ 
tion of our halos is lower than what is often assumed in the liter¬ 
ature. When baryon effects are neglected, an inner slope close to 
7 = 1.26 is clearly incompatible with the simulations. 


3.2 Profiles in the simulations with baryons 

We now turn towards the dark matter profiles in the simulations 
including baryons. In Fig. 2, we show the dark matter density 
profiles of the four halos simulated with the full baryonic model. 
As for the previous figure, the lines are thin at radii less than 
the convergence radius of the simulation Rpoa and the dashed 
lines correspond to the NFW and gNFW profiles normalised to 
p(r©) = 0.4 GeV • cm“®. The simulated profiles present two in¬ 
teresting features when compared to the DMO results. In the range 


Table 1. Properties of the four simulated DMO halos (Fig. 1) and the best¬ 
fitting NFW parameters to their dark matter density profiles. The spherical 
overdensity masses and radii are given with respect to the critical density of 
the universe. 


Halo 

M 200 

[Mq] 

R 200 

[kpc] 

Rpos 

[pc] 

[kpc] 

PDM(r©) 
[GeV • cm“' 

1 

1.65 X IOI 2 

243.2 

435 

22.4 

0.290 

2 

1.09 X 10^2 

212.0 

445 

20.1 

0.132 

3 

1.35 X 10^2 

226.9 

344 

23.2 

0.162 

4 

1.39 X I 012 

229.4 

358 

19.8 

0.281 



Figure 2. The dark matter density profiles of the four halos in the simula¬ 
tions with baryons physics (red solid lines). Thinner lines ai'e used at radii 
smaller than the convergence radius of the simulation. The vertical dotted 
line indicates the simulation’s gravitational softening length. The green dot- 
dashed and blue dashed lines correspond to an NFW and gNFW with 7 = 
1.26 profiles respectively, bolh normalised to p(rQ) = 0.4 GeV ■ cm“® 
and with a scalelength Tb = 20 kpc. The profiles display a logarithmic 
slope steeper than —1 between r Ki 1.5 kpc and r Ri 10 kpc, in broad 
agreement with the profiles inferred from observations by Galore et al. 
(2015b). At radii r 1 kpc the profile is significantly shallower than NFW 
profiles are. 

1.5—10 kpc, the profile is significantly steeper than an NFW profile 
and at radii r < 1 kpc, the profiles display a significant flattening. 
Our profiles thus display a combination of dark matter contraction 
and a flattening further in. The properties of the dark matter distri¬ 
bution in the central regions are, however, particularly sensitive to 
the choice of subgrid model so these results, particularly the flatten¬ 
ing of the profile, should be regarded as tentative and, by no means, 
as a generic prediction of ACDM. 

These halos are clearly not well described over their entire ra¬ 
dial range by any of the profiles commonly found in the literature. 
The main properties of the halos are given in Table 2. Note that 
in agreement with the findings of Schaller et al. (2015a) for the 
lower-resolution periodic EAGLE volume, the halo masses M 200 
(and hence radii i? 20 o) are lower than in the simulation that did 
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Table 2. Properties of the four simulated EAGLE halos (Fig. 2) and the best¬ 
fitting gNFW asymptotic slope 7 to their dark matter density profiles in the 
radial range r > 1.5 kpc. 


Halo 

AT 200 

[Mq] 

R 200 

[kpc] 

Rp03 

[pc] 

7 

[-] 

PDM(rQ) 
[GeV ■ cm“' 

1 

1.56 X 10^2 

238.8 

559 

1.38 

0.310 

2 

1.01 X 10^2 

206.8 

592 

1.47 

0.160 

3 

1.12 X 10^2 

213.7 

438 

1.73 

0.204 

4 

1.16 X I 012 

216.2 

462 

1.49 

0.280 


not include baryon physics. A consequence of the steepening of the 
profile due to contraction is the slight increase in the local dark mat¬ 
ter density pDM(f 0 ) (column 6 of the tables), which, however, re¬ 
mains lower than the commonly adopted value of 0.4 GeV • cm“® 
in each case. Clearly, the simulated profiles will not be well de¬ 
scribed by a gNFW profile at radii r < 1.5 kpc. It is, however, in¬ 
structive to find the best-fitting profile at larger radii for comparison 
with the models used to characterise the Fermi excess. The best¬ 
fitting asymptotic slopes are given in column 5 of table 2. For all 
four halos, the slopes are steeper than the value (7 = 1.26 ± 0.15) 
found by Galore et al. (2015b) when modelling the GC excess. We 
note, however, that the simulated profiles are in broad agreement 
with the gNFW profile of Galore et al. (2015b) (blue dashed line), 
if the overall normalisation is, once again, ignored. The baryons 
significantly steepen the profiles at radii r > 1.5 kpc. 

At radii r < 1.5 kpc, the density profiles deviate significantly 
from the cusp seen in the DMO simulation. At the resolution limit, 
Rpq 3 ~ 450 — 600 pc, the simulated profiles exhibit a density 
between 2.5 and 4.2 times lower than the best-fit gNFW profile 
inferred from observations. This flattening is an important feature 
since the densest regions of the haloes dominate the 7 -ray emission. 
No such flattening was seen in the reference EAGLE simulations, 
which, however, had over 200 times poorer mass resolution than 
the Local Group APOSTLE simulations used here. Here, the flatten¬ 
ing is well resolved since it occurs at radii significantly larger than 
the Power et al. (2003) convergence radius. This indicates that the 
flattening is not a result of poor sampling but rather a consequence 
of our specific choice of subgrid model. 

At high redshift all four examples had developed the cuspy 
central profile that is characteristic of dark matter haloes. How¬ 
ever, the inner profiles flattened during events that are clearly as¬ 
sociated with violent star formation in the inner few kiloparsecs. 
In one example, the cusp regenerated before being flattened again 
by a new episode of violent star formation activity. This phe¬ 
nomenon is reminiscent of the cusp-destroying mechanism pro¬ 
posed by Navarro et al. (1996a) in which the sudden removal of 
dense, self-gravitating gas from the centre by a starburst redis¬ 
tributes the binding energy of the central regions. Related processes 
have been seen in recent simulations, mostly of dwarf galaxies 
(e.g. Govemato et al. 2010; Duffy et al. 2010; Maccio et al. 2012; 
Pontzen & Govemato 2012; Teyssier et al. 2013; Onorbe et al. 
2015). We defer a detailed discussion of the causes of this inter¬ 
esting phenomenon to a separate study. 

3.3 Resolution and convergence considerations 

In the previous two subsections, we used the criterion of Power 
et al. (2003) as the resolution limit of our simulations. This crite¬ 
rion, based on the two-body relaxation timescale, was derived using 
purely collisionless simulations and was designed to ensure that 



Figure 3. Convergence test for the dark matter profiles of one of the haloes 
(Halo 1) in both the simulation with baryons (upper set of lines) and with¬ 
out baryons (lower set of lines, rescaled by a factor of 0.2 for clarity). The 
green, blue and red lines correspond to the simulations with a dark matter 
particle mass of 7 X 10® Mq, 6 X 10® Mq and 5 X lO'* Mq respec¬ 
tively. Dashed lines are used at radii smaller than the convergence radius, 
Rp 03 . of each resolution. The vertical dotted lines indicate the softening 
length of each simulation. The thick yellow lines show an NFW profile 
with Tb = 20 kpc. The profiles in the DMO simulation are well converged 
even at r < iJpo 3 . In the EAGLE simulation, the profiles display a lower 
level of convergence but nevertheless agree within 30% at r > Rpo 3 . The 
differences between the various resolutions are, however, smaller than the 
difference relative to the DMO case. The behaviour of this particular halo is 
typical of the four cases we simulated. These are all shown in Appendix A. 

the enclosed mass at a given radius, i?po 3 , is within 10 % of the 
value obtained using a higher resolution simulation. In the cases 
where baryonic effects are simulated, it is unclear whether this cri¬ 
terion still applies and even whether numerical convergence in the 
usual sense can be achieved (see discussion in Schaye et al. (2015)). 
Schaller et al. (2015a) demonstrated that stacked haloes, extracted 
from the EAGLE volumes, are well converged when using this sim¬ 
ple criterion but it is unclear whether this remains tme when in¬ 
dividual haloes are considered. Of particular concern is the use of 
the pressure floor for the densest gas, which sets an artificial scale 
below which gas cannot be compressed. For our simulations, at 
all resolutions, this pressure floor ensures that Jeans lengths above 
Aj « 750 pc can be resolved and prevents the collapse of gas 
clouds of smaller sizes. It is therefore possible that the profiles may 
be modified by this pressure floor at radii r ~ Aj, which, inci¬ 
dentally, is similar to the value of i?po 3 in our highest resolution 
simulation. 

In Fig. 3 we show the density profiles (multiplied by to 
reduce the dynamic range) of the dark matter extracted from one 
of our simulations run at three different resolutions, separated by 
factors of 12 in particle mass. The bottom set of curves are ex¬ 
tracted from the simulation without baryons (and have been mul¬ 
tiplied by a factor 0.2 for clarity), whilst the upper set of lines are 
taken from the simulations with baryons. Dashed lines are used at 
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radii r < Rpoa and the vertical dotted lines indicate the softening 
lengths in each of the three simulations. As a guide, NFW profiles 
with ra = 20 kpc and the same value of pDM(rQ) as the highest 
resolution simulation are shown using yellow lines. Similar figures 
for the three other haloes are shown in Appendix A. 

In the DMO simulations (lower set of curves), the profiles of 
different resolution converge very well and, as expected, the cri¬ 
terion of Power et al. (2003) (which refers to the enclosed mass) 
provides a good, conservative, estimate of the radius above which 
the density profiles are converged. In fact, the density profiles are 
converged to r > 0.7i?po3. The profiles are very well described 
by the NFW functional form and show an asymptotic inner slope 
of —1. The deviation from NFW seen at r 30 kpc is due to the 
presence of substructures in the haloes. 

The situation is different for the haloes with baryonic physics 
(upper set of curves). At r < 10 kpc, the profiles show significant 
differences when the resolution is varied. Important differences are 
especially visible between the two highest resolutions (blue and 
red curves). Clearly, the criterion of Power et al. (2003) is no longer 
suitable because the profiles show differences even at r > J?po 3 but 
note that in that range all three resolutions nevertheless agree within 
30%. Despite this poorer level of convergence, the profiles display 
similar general trends across resolution levels. The profiles are sig¬ 
nificantly steeper than NFW at 1.5 kpc < r < 10 kpc and sig¬ 
nificantly shallower than NFW at smaller radii. These differences 
relative to the NFW profile are larger than the differences between 
the various resolution simulations, suggesting that the trends seen 
are a generic outcome of the subgrid model assumed in our simu¬ 
lation even if their exact magnitude is difficult to establish. In the 
remainder of this paper, we will restrict our analysis to the high¬ 
est resolution simulations, but these limitations should be borne in 
mind when interpreting our results and evaluating the generality of 
our conclusions. 


3.4 Sphericity of the distribution 

In order to characterise the morphology of the dark matter annihila¬ 
tion signal at the centre of the Milky Way, it is interesting to study 
the shape of the dark matter distribution. The profiles described so 
far assumed a spherically symmetric dark matter density profile. 
With the higher precision of the measurements of the excess and 
the increasing understanding of the GDE, it will soon be possible 
to measure deviations from a perfect sphere. For instance, the pres¬ 
ence of a “dark disc” (Read et al. 2008) would enhance the signal 
in the plane of the galactic disk and hence break the symmetry of 
the signal. This would also make the signal more difficult to dis¬ 
entangle from astrophysical components associated with the disk, 
such as point sources. 

In order to test this, we computed the inertia tensor of our 
four halos using all the dark matter within a spherical aperture of 
1 kpc from the centre. This distance corresponds to « 7° on the 
sky and hence encompasses the majority of the 7 -ray flux in the 
direction towards the Galactic Centre that would result from dark 
matter annihilation in the MW. We then compute the three eigen¬ 
values a > 6 > c of the inertia tensor and report the values in Table 
3 for both simulations with and without baryons. 

As can be seen, the axis ratios are very close to unity, indi¬ 
cating only very small deviations from sphericity and hence no ob¬ 
vious anisotropy feature in the signal. We also find no alignment 
between the main axis of the dark matter distribution in the inner 
1 kpc and the plane of rotation of the stars. It is interesting to note 
that the simulation with baryons yields more spherical distributions 


Table 3. Axis ratios (a > b > c) infeiTed from the inertia tensor of the 
matter within 1 kpc from the centre of the galaxies for both the halos in the 
DMO and EAGLE local group simulations. 

DMO EAGLE 


Halo 

hja 

cjh 

bfa 

cjh 

1 

0.888 

0.973 

0.989 

0.953 

2 

0.863 

0.957 

0.975 

0.971 

3 

0.850 

0.984 

0.981 

0.941 

4 

0.879 

0.964 

0.987 

0.988 


close to the centre than its counterpart without baryons. We verified 
that repeating the exercise with apertures of 0.5, 2 and 3 kpc yields 
similar results. 


3.5 Position of the centre 

Another potential source of systematics in the analysis of the GC 
excess is the position of the centre of the dark matter distribu¬ 
tion. If the highest-density part of the dark matter profile is offset 
from the centre of the stellar distribution, then this offset should be 
detectable in observations. In their simulation of a single Milky 
Way-like halo, Kuhlen et al. (2013) found a sizeable offset of 
300 — 400 pc between the centre of the stellar distribution and 
the peak of their dark matter distribution. If such an offset was in¬ 
deed present in the Milky Way, then an offset of ~ 2° between 
the GC and the peak of the dark matter annihilation signal should 
be visible. In their study based on the EAGLE simulations, Schaller 
et al. (2015b) found that the offset between the peak of the dark 
matter density distribution and the centre of the stellar distribu¬ 
tion is typically smaller than the softening length of the simulation 
(e = 700 pc in their case). Repeating their analysis on our four 
simulated Milky Way halos, we find offsets between 22 and 43 pc, 
well below the size of the softening length (e = 134 pc), indicating 
that the offsets are consistent with zero. For all practical purposes 
and given the current resolution of instruments, the centre of the 
dark matter distribution is hence coincident with the centre of the 
stellar distribution. 


4 DARK MATTER ANNIHILATION SIGNAL 

Now that the dark matter profiles have been characterised, we turn 
to the derivation of the corresponding annihilation signal as ob¬ 
served by a virtual instrument located at the position of the Sun 
and pointing towards the centre of the Milky Way. 


4.1 J-factor for the simulated halos 

In the case of a dark matter particle that annihilates into photons or 
into particles that generate photons in their decays, the photon flux 
(in units of GeV“^ • cm“^ • • sr“^) at a given angle, 4^, on the 

sky away from the GC is given by 


dN 

'dE 


(vF) 


{av) dN-y 
8nm'^ dE 




( 2 ) 


where is the mass of the dark matter particle, (av) is its ve¬ 
locity averaged total annihilation cross section, dN-y/dE is the av¬ 
eraged energy spectrum of photons produced per annihilation and 
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/(^') is the integral along the line of sight of the square of the dark 
matter density. This so-called “J-factor” reads 

/(«-)=/ p!,M(r(s,^))ds, (3) 

J 1.0.S. 

with the variable s running along the line of sight axis from s = 0 
to s = oo and 

r(s,T') = ^(rQ — s cos Vl')^ -|- (s sin (4) 

giving the distance from the GC for a particular angle on the sky 
T' and distance to the GC, tq. The differential intensity dN/dE is 
hence the product of the J-factor, given by the distribution of dark 
matter, and the particle physics model assumed. As a consequence, 
within a reasonable range, the precise normalisation of the J-factor 
is irrelevant since a similar signal can be recovered by altering the 
particle physics model. To simplify the comparison with the anal¬ 
ysis of the GC excess, we have, thus, normalised our simulated 
profiles such that pDM(r' 0 ) = 0.4 GeV • cm~^. 

In the top panel of Fig. 4 we show the J-factor (Eq. 3) as a 
function of galactic latitude b (at galactic longitude / = 0 °) for our 
four simulated profiles, normalised to the same local dark matter 
density. The red and yellow lines correspond to the the dark matter 
profiles in the simulations with and without baryons respectively. 
The green dot-dashed and blue dashed lines correspond to NEW 
and gNEW with 7 = 1.26 profiles respectively with a scale length 
Tg = 20 kpc and the same normalisation local dark matter den¬ 
sity as our normalised halos. The shape of the J-factor profile is 
different in the runs with and without baryons. The contraction of 
the dark matter due to baryons increases the J-factor by a factor of 
~ 2 at angles fe > 4° from the GC, when compared to an NEW 
halo. In that angular range, the J-factor is also larger than the one 
obtained for a gNEW with a slope 7 = 1.26 (Calore et al. 2015b). 
Closer to the GC, the simulated J-factors display a shallower slope 
and values lower than the gNEW model. 


4.2 Extrapolation of the profiles towards the centre 

As most of the dark matter annihilation signal originates from the 
inner few hectoparsecs, it is necessary to extrapolate our findings 
from section 3.2 to smaller radii. As the annihilation signal in¬ 
creases with the square of the local density, one can ask what the 
highest density is that can be reached in the inner regions given 
the constraints on the density and enclosed mass at the smallest 
converged radius. Assuming that the profile is not hollow and that 
the logarithmic slope is monotonic, it is straightforward to show 
that the only asymptotic power law that can be used to extrapolate 
the profiles from a given radius r towards the centre has a slope 
7max = 3(1 — 47 tr®p(r)/ 3 M(r)) (Navarro et al. 2010). Setting 
r to the convergence radius of the halos, i?po 3 , we obtain slopes 
in the range 7 max = 0.55 — 1.22 for our four halos. The J-factors 
resulting from these extrapolations are shown on Fig. 4 using thin 
lines. They allow us to set upper bounds on the J-factor for angles 
6 < 3°. Even with this power-law extrapolation, the flux is lower 
than the gNEW profile with a slope 7 = 1.26. 


4.3 Gamma-ray flux morphology 

In the bottom panel of Fig. 4, we show the emission at = 2 GeV 
for our J-factors, assuming the best-fitting particle physics model 


of Calore et al. (2015b)®. For comparison, we show the flux in¬ 
ferred from the GC excess by Hooper & Goodenough (2011),Bo- 
yarsky et al. (2011), Gordon & Macias (2013), Hooper & Slatyer 
(2013), Daylan et al. (2014) and Abazajian et al. (2014) with the 
error bars indicating the ilcr statistical error (not shown when 
smaller than the symbols). The observed intensities were rescaled 
following the procedure highlighted in Calore et al. (2015c), tak¬ 
ing into account the assumed excess profiles. Note that individual 
measurements are more than 3cr discrepant with each other. The 
green shaded regions indicate the best-fitting model of Calore et al. 
(2015b). Their model assumed a gNFW profile for the dark matter 
profile and used 60 GDE templates in their likelihood analysis of 
10 regions of interest on the sky located around the galactic centre. 
The width of the green regions on the figure indicates both the sta¬ 
tistical uncertainty and the posterior range of the GDE modelling 
around the best-fitting gNEW profile. The uncertainty on the slope 
of the profile is not shown. A similar analysis, performed by Calore 
et al. (2015c), of the preliminary results of the Fermi collaboration 
is shown as a yellow shaded region. The grey shaded region on the 
left of the plot indicates the radial range over which the emission 
from the ISM gas dominates the GDE models (Calore et al. 2015c). 
Similarly, the grey shaded region at the bottom of the panel indi¬ 
cates the level of 7 -ray flux expected from the extended “Fermi 
Bubbles” (Su et al. 2010), thought to be the remnant of past AGN 
activity. We use the extrapolation, assuming a constant density, to 
lower latitudes of the flux estimated by Ackermann et al. (2014). 
The flux originating from the annihilation of dark matter is higher 
than the contribution of the Fermi Bubbles at angles b < 15°, mak¬ 
ing the radial range 2° < b < 15° ideal for the study of the excess 
(Calore et al. 2015c). The resolution of our simulations is, hence, 
well matched to this requirement. 

Our simulated profiles (red lines) are in good agreement with 
the 7 -ray data for angles & > 3°. This is expected since over the 
relevant radial range, the profiles have a similar form to the gNFW 
profile with asymptotic inner slope, 7 = 1.2 — 1.3, inferred di¬ 
rectly from the data (assuming that the emission is due to dark mat¬ 
ter annihilation). At smaller angles, three of the four extrapolated 
profiles give significantly less emission than observed, whilst the 
fourth is in good agreement with the data. We stress, however, that 
this extrapolation gives the largest power-law signal at the GC. The 
predicted emission at fe < 2 ° could be boosted by adjusting the 
particle physics model but there is a danger that such adjustments 
could lead to an overprediction of the emission at larger angles. 


5 SUMMARY & CONCLUSION 

In this study we investigated the dark matter density profiles of 
four Milky Way galaxies simulated using a state-of-the-art hydro¬ 
dynamics code and subgrid model. We specifically focused on the 
inner few kiloparsecs of the dark matter distribution in order to 
refine earlier predictions for the properties of dark matter annihila¬ 
tion emission. The careful treatment of baryons in our simulations 
allows us to understand and analyse the effects that baryons have 
on the dark matter distribution. These are not negligible and give 
rise to haloes whose properties differ significantly from those of 
the “standard halo model” often used in dark matter detection stud¬ 
ies. Whilst our simulations are among the highest resolution exam¬ 
ples of their kind currently available, we are only able to establish 

® = 46.6 GeV, (av) = 1.60 X 10“^® cm® for a bb final 

annihilation state. 
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Figure 4. Top panel: The J-factor as a function of galactic latitude b inferred from our four simulated halos both with (red lines) and without (yellow 
lines) baryon effects. The green dot-dashed and blue dashed lines correspond to NFW and gNFW with 7 = 1.26 profiles respectively with a scale length 
rs = 20 kpc. To ease comparison, all profiles have been normalised to yield PDM(rQ) = 0.4 GeV ■ cm“®. The thin lines correspond to power-law 
extrapolations of our simulated profiles (see text). The scale at the top indicates the minimum radius intersected by a line of sight at the given galactic latitude 
b. Bottom panel: Emission at i? = 2 GeV for our halos assuming the best-fitting particle physics model from Galore et al. (2015c). Data points with error 
bars show the best-fitting models of Hooper & Goodenough (2011), Boyarsky et al. (201 1), Gordon & Mactas (2013), Hooper & Slatyer (2013), Daylan et al. 
(2014) and Abazajian et al. (2014) to the Fermi GeV excess. The green shaded regions indicate the excess emission and its statistical uncertainty for a fixed 
gNFW profile derived by Galore et al. (2015b) and the yellow shaded region indicates the preliminary results of the Fermi-LAT team. The vertical grey shaded 
region indicates the radial range where uncertainties in the GDE modelling due to tt® emission from HI and H 2 regions dominate the foreground templates 
used in the analysis of the data (adapted from Galore et al. 2015c). Similarly, the shaded region at the bottom indicates the flux intensity of the Fermi bubbles. 


convergence within a tolerance of 30% across different resolutions 
over the radial range of interest (sec. 3.3). We feel that this level 
of convergence is sufficient to support our conclusions but higher 
resolution simulations will be needed to test this supposition. 

We can summarise our findings as follows: 

• As seen in previous simulations (e.g. Dubinski 1994; Abadi 
et al. 2010; Bryan et al. 2013), the central concentration of baryons 
significantly reduces the asphericity typical of haloes in dark 
matter-only simulations. The distribution of dark matter in the 
inner 500 pc is very close to spherical with axis ratios, h/a > 0.96, 
in all cases. 

• There is no detectable offset between the position of the GC and 
the peak of the dark matter distribution. The largest offset found in 
our halos is 45 pc, much smaller than the softening length of the 


simulations (e = 134 pc). 

• The condensation of baryons at the halo centre causes the halo 
to contract slightly. The halo density profiles end up having steeper 
profiles than NFW in the radial range r = 2 — 10 kpc. 

• In the inner 1.5 kpc, which are well resolved in our simulations 
the dark matter halo density profiles develop significant flattening. 
This feature is likely to be associated with violent star formation 
events that take place during the early stages of galaxy formation. 
It must be borne in mind that effects of this kind are sensitive to 
the specific subgrid physics model and, at this point, they must not 
be regarded as a generic prediction of the ACDM model. 

• The predicted dark matter annihilation emission signal is in good 
agreement with the detection of extended 7 -ray emission in excess 
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of the known foregrounds by the Fermi satellite at galactic latitudes, 
6 > 3°, where our haloes are well resolved. A simple extrapolation 
of the density profiles in our simulations to smaller angles predicts 
a 7 -ray flux significantly lower than is measured, in three of the 
four cases, suggesting possible contributions from other sources to 
the excess. In the fourth case, the annihilation signal from the ex¬ 
trapolation is in broad agreement with the reported measurements 
close to the Galactic Centre. 

The analysis of the Fermi excess has so far been performed 
assuming a gNFW profile or other parametric profile forms for the 
dark matter. Future, more precise studies, would benefit from using 
the more realistic profile shapes derived directly from hydrodynam- 
ical simulations. This should help disentangle the signal from dark 
matter annihilation from the galactic diffuse emission. 
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Figure Al. Same as Fig. 3 but for all four halos considered in this study. The DM density profiles in all four haloes agree to better than 30% at r < 10 kpc. 
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